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Abstract 

The MACRO experiment at Gran Sasso provides means for detailed studies of 
multiple coincident penetrating cosmic ray muons. In this paper we concentrate 
on the studies of the ultrahigh energy primary cosmic ray composition using muon 
bundle multiplicities, muon pair lateral and angular separation distributions. 



1. INTRODUCTION 

The chemical composition and energy spectra of primary cosmic rays bombarding the 
top of the Earth's atmosphere contain information about the origin of cosmic rays, pro- 
viding constraints on the modeling of acceleration and propagation in the interstellar 
medium. Experiments with balloons and satellites provide direct measurements of the 
composition and energy spectra of the primary species up to energies of several hundred 
TeV. At higher energies, direct measurements require impractically large detectors and/or 
long exposure time due to the steeply falling energy spectra. Consequently, very little is 
known about the composition of primary cosmic rays at energies greater than several 
hundred TeV [1]. However, indirect information on their composition can be extracted 
from experiments with secondary particles from atmospheric cascades, such as surface 
measurements of extensive air showers or underground studies of penetrating high energy 
muons. 

The rates of underground muon bundles of different multiplicities, as well as their lateral 
and angular separations, are sensitive to the chemical composition and energy spectra of 
the primary cosmic rays, above a threshold determined by the rock overburden (~ 50 TeV 
at Gran Sasso). This sensitivity arises from the fact that heavy nuclei tend to generate 
a larger yield of charged pions and kaons, in the forward fragmentation (high rapidity) 
region, than that of light nuclei. It is these pions and kaons which can decay into the 
penetrating muons observed deep underground. 

In addition to the composition of cosmic rays, the properties of hadronic interactions 
at very high energies and in the high rapidity region will also affect the multiplicity and 
separation of muons observed underground. Interactions with higher overall cross sections 
will tend to cause primaries to interact higher in the atmosphere which can result in wider 
spread between secondary muons. Similarly, heavy nuclei will tend to interact higher in 
the atmosphere and will also result in larger muon spread. However, if an interaction 
tends to produce high multiplicity with little increase in transverse momentum, then 
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multiplicity will rise if the energy is sufficiently high (it can fall below a particular energy 
threshold) but there will be no increase in transverse momentum. Hence, there are a 
number of ways that the primary energy, composition, and interactions can influence 
the muon distributions. By studying these distributions and their correlations, various 
composition and interaction models can be constrained. 

In this paper we report results from analyses [2] of the multiple coincident penetrating 
cosmic ray muon data collected by the MACRO detector. We have measured the muon 
bundle rate as a function of multiplicity as well as the muon pair lateral and angular 
separation distributions. These measured quantities are compared with Monte Carlo 
predictions for different primary cosmic ray composition models and interaction models. 



2. THE MACRO DETECTOR 

The MACRO (Monopole, Astrophysics, and Cosmic Ray Observatory) detector [3], a 
large area, deep underground detector, was fully completed recently at the Gran Sasso Na- 
tional Laboratories in central Italy. It consists of large liquid scintillator counters, limited 
streamer tubes, and plastic track-etch detectors, which offer three redundant techniques 
for its primary physics goals of monopole searches [4]. Its large size (76 X 12 X 9 m"^) 
and excellent tracking capability allow detailed studies of multiple coincident penetrating 
cosmic ray muons. The rock overburden has a minimum thickness of 3200 meters of water 
equivalent, setting the surface muon energy threshold of ~ 1.4 TeV. Above this energy, 
the probability for a muon to penetrate this minimum thickness is 30% or greater. 

The full MACRO detector including its upper deck started to take muon data in March 
1994. This paper uses data collected from the lower deck, which commenced data taking 
in June I99I. Consisting of six supermodules, the lower deck measures 76 X 12 X 4.8 m"^. 
It is surrounded on all sides by planes of large liquid scintillator counters. The tracking 
system consists of ten horizontal layers of streamer tubes separated by 32 cm of crushed 
rock absorber. Each tube is 12 m long, 3x3 cm^ in cross section, and utilizes a 100 /im 
anode wire and a graphite cathode. A two-dimensional readout is performed using signals 
from anode wires and external 26.5° stereo angle strips. Six additional vertical layers of 
streamer tubes cover each vertical side. The intrinsic angular resolution is 0.2° for muons 
crossing ten horizontal planes. Taking into account multiple coulomb scattering in the 
rock overburden, the overall angular resolution is estimated to be ~ 1°. 

MACRO accumulates underground muon data at the rate of ~ 6.6 X 10^ events/live 
year. Approximately ~ 4.0 X 10^ events/live year exhibit multiple muon tracks with 
lateral separations up to 70 m and ~ 1.6 X 10"^ events/live year have multiplicities of ten 
or greater. Combined with our ever improving ability to model cosmic ray showers, these 
data provides a unique opportunity to study ultrahigh energy cosmic rays. 
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3. COSMIC RAY COMPOSITION 

We have studied the primary cosmic ray composition using the multi-muon event rate 
as a function of multiphcity. This analysis uses a data sample collected in 3295 hr live 
time with all six lower supermodules. This data sample contains ~ 2.5 X 10^ muon events, 
of which ~ 1.5 X 10^ are multiple muons. The event selection uses the criteria established 
in a previous analysis [5], which used data from only two supermodules. In the previous 
analysis, a considerable amount of visual scanning was performed to establish the actual 
multiplicity. In this analysis, we correct the reconstructed multiplicity using a GEANT- 
based [6] detector simulation program. The following physics and detector effects are 
taken into account: electromagnetic showering down to 500 keV, charge induction of 
the streamer signal onto the stereo strip, electronic noise, inefficiencies, and failures of 
the tracking algorithm especially for high multiplicity events, track shadowing at small 
separations, etc. The simulated data are used to calculate the correction factors for 
transforming the reconstructed multiplicities in the two projective readout views into 
an actual multiplicity. This allows an objective assignment of the muon multiplicity, 
reducing the systematic uncertainties in the previous analysis dominated by the scanning 
uncertainties. Figure 1 shows the muon bundle rates for the one, two, and six supermodule 
data samples. The increase of acceptance with detector size is reflected in the increase in 
muon rates and sampling of high multiplicity events. 

The experimental data are compared with the results of full simulations of the primary 
interaction, the atmospheric cascades, the muon propagation in the rock, and the afore- 
mentioned detector response. Two different shower simulation codes are used: HEMAS 
[7] with the addition of nuclear fragmentation based on the semi-superposition model [8] 
and SIBYLL [9]. The simulations are described in more detail below in Section 4. This 
paper reports results based on the HEMAS simulations. 

We considered three different primary composition models: light and heavy composi- 
tions used in our previous analysis [5], and a constant mass composition (CMC) with fixed 
spectral indices [10]. The light and heavy compositions are extreme models: at increas- 
ingly higher energies the light composition contains a large proton component while the 
heavy composition contains a large Ee component. However, the models are constrained 
to reproduce the known abundances and spectra directly measured at < 100 TeV and to 
agree with extensive air shower measurements at higher energies. Therefore a comparison 
of the muon experimental rates with the predictions of these models gives an indication 
of the sensitivity of the experiment to primary composition. 

Eigure 2 shows the calculated energy ranges of the primary cosmic rays at a detection 
efficiency of 90% as functions of detected multiplicities for light and heavy composition 
models. The plots demonstrate that the primary energies explored this way increase with 
muon multiplicity. In particular, events with detected multiplicity >, 10 originate from 
primaries in an energy region entirely above the "knee" , the steepening point of the cosmic 
ray energy spectra at ~ 100 TeV. 
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Figure 1. Muon bundle 
rates as functions of muon 
multiplicity for the one, 
two, and six supermodule 
data samples. 
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Figure 2. Calculated en- 
ergy ranges of the primary 
cosmic rays at a detection 
efficiency of 90% as func- 
tions of detected multiplic- 
ities. The bold curves indi- 
cate the mean primary en- 
ergies. 
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Figure 3 shows the dependence of the average primary mass (A) on energy for the 
three composition models described above, as well as for the SIGMA model which is 
based on fits to direct measurements from 1-100 TeV and extrapolated up to the "knee" 
region [11]. Therefore the SIGMA model can be viewed as a reference (A), following the 
energy dependence of direct measurements at low energy. As one can see from this figure, 
the light and CMC models are in the same range of (A) as the direct measurements. The 
heavy model, on the other hand, has a completely different energy dependence in nearly 
the entire energy range relevant to MACRO multi-muon events. 

The ratios between the rates predicted by the HEMAS simulations and the experimental 
data are shown in Figure 4. One feature of this figure is that the measured multi-muon 
rates at low multiplicities [N^ < 4) are higher than those predicted by the Monte Carlo 
regardless of composition model. These events originate in primaries with energies less 
than a few hundred TeV, for which the three models are very similar and in agreement 
with direct measurements. We investigated the effects of our present uncertainties in 
the rock overburden and muon propagation as possible sources of this disagreement, and 
found from Monte Carlo simulations that while these uncertainties affect the absolute 
muon rates they do not affect the shape of the multiplicity distribution. An analysis 
based on this shape is in progress. 

Figure 4 shows that our data are inconsistent with the predictions of the heavy compo- 
sition at high multiplicity and favor a lighter model. Therefore our data do not favor the 
hypothesis used in the heavy model which requires a dramatic change of primary compo- 
sition toward the pure iron element immediately above the "knee." The data provide a 
better fit to models with fiat or slowly increasing (A) as a function of primary energy, as 
in the light and CMC models. 



4. MUON DECOHERENCE 

The lateral separation of underground muon pairs has been demonstrated to be sensitive 
to hadronic interaction models, as well as to the primary composition models, allowing 
the rejection of some simplified cascade treatments [12]. Here we report an analysis of a 
data sample of ~ 5.8 X 10^ muons and ~ 7600 hr live time, in which 1.9 X 10^ muon pairs 
are reconstructed using the criteria described in [13]. 

The muon pair lateral separations are traditionally analyzed using a detector- 
independent "decoherence function" [14] defined as the rate of muon pairs per unit area, 
per steradian, per pair separation determined on a plane orthogonal to the pair direction: 

1 fd'N,{r,e,cf>)/drdn^^ 



A{e,<p) 



where Np/ drdil is the density of pairs at distance r and incidence angle (^, (/>), A is the 
projected detector area in the (^, (/>) direction, is the total solid angle defined by the 
limits of integration, and T is the exposure time. In this analysis, the two independent 
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methods described in [12] are used to compute the decoherence function and they yield 
the same results. 

To compare with composition and interaction models, a detailed Monte Carlo cal- 
culation is required to simulate the production and propagation of muons through the 
atmosphere and mountain overburden. The previous work [12] used the parameterized 
results of the HEMAS code [7] to generate cosmic ray showers, in order to save computer 
time. The parameterized formulae gives the number and the spatial distribution of the 
underground muons as functions of primary mass, energy, and direction, but it ignores 
the correlation between muon multiplicity and lateral distribution. Therefore we have 
chosen in this new analysis to perform full simulations of the atmospheric cascade and 
muon propagation. In preparing the full HEMAS based Monte Carlo, we investigated 
several possible systematic effects on the simulation, including: (1) biases due to the fail- 
ure of the parameterization to account for the correlation between muon multiplicity and 
lateral distribution; (2) uncertainties in the primary interaction cross section; (3) different 
possible models of nucleus-nucleus interactions (superposition vs. fragmentation); (4) the 
treatment of energy loss and multiple scattering in the rock; and (5) the effect of the 
geomagnetic field on cascade development. 

The results generally showed that the dependence of the underground muon separation 
on the details of the Monte Carlo was weak. A comparison of the full HEMAS based 
Monte Carlo to the parameterized version using the constant mass composition [10], for 
instance, showed that the average muon pair separation increases by only 2% in the 
more complete version. A 10% increase in the primary interaction cross section raises 
the muon production height by 3%, and consequently increase the average underground 
muon separation. There is a great deal of uncertainty in the nucleus-nucleus interaction 
model, due to the lack of accelerator experimental data for the energies of interest. The 
HEMAS code handles cascades generated by heavy nuclei (mass A, total energy E) in 
the superposition scheme, as A independent nucleons of energy E/A. Replacing this 
with a more realistic model [8] causes larger fiuctuations in shower development, but the 
decoherence distribution is not affected within the present statistics of the simulated data. 
HEMAS muon transport through the rock was compared to that implemented in GEANT 
[6], a more recent simulation developed to model high energy accelerator events; we found 
no noticeable difference as far as the lateral spread of muons is concerned. Einally, an 
increase in the average muon separation of about 5% is achieved by considering the effect 
of the geomagnetic field on shower development. 

In Eigure 5 the experimental pair separation distribution for the entire data sample is 
compared with Monte Carlo predictions for two extreme (light and heavy) composition 
models. In the previous analysis [12] of data from only two supermodules, we found a good 
agreement with the Monte Carlo predictions up to the maximum attainable separation 
(~ 20 m). With all six lower supermodules the measured decoherence function at large 
separations is higher than that from the simulations. The average separation is 10.9ib0.2m 
for the real data, 10.5 m for the heavy model, and 9.4 m for the light model. The measured 
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Figure 5. The decoherence 
function function of the 
muon pair separation. 



separation is larger than the one from the simulations. We have also investigated the 
dependence of muon pair separation on rock depth and zenith angle. Tables 1 and 2 
summarize the main results. It is worth noting that the experimental and simulated values 
follow the same behavior and are close to one another. But the experimental values are 
larger than the simulated ones, following the same trend in the overall distribution. 

Since the heavy model is a disfavored extreme model (see Section 3), Figure 5 may indi- 
cate that the shower development is not yet treated properly in the Monte Carlo and that 
the hadronic interaction model presently used for this analysis may need further improve- 
ment. In particular, we are investigating factors which affect the transverse momentum 
distribution, including possible nuclear effects. 

In order to better understand the role of the hadronic interaction model, we have 
performed a comparative simulation with the SIBYLL code [9]. No substantial changes, 
however, were observed as far as muon separation is concerned. Such a model in fact 
predicts a slightly lower average muon separation (on the order of 10% less). In the 
future, other models will be considered, such as the DPMJFT code [15], which has a 
more complete treatment of nuclear effects than does SIBYLL. 

5. MUON DECORRELATION 

Detailed measurements of quantities relative to underground multiple muon events can 
provide information on primary interactions and muon propagation through the rock over- 
burden. In particular, the differential distribution of spatial and angular separations [x 
and (j)) in multiple muon events, dN j dxdcf), is sensitive to the physics of muon production 
and propagation. It provides information on the total primary cross section as well as to 
the transverse momentum distribution of the parent hadrons of the underground muons. 
It also provides a measure of the effects of muon interactions in the rock overburden. 



9 



Table 1 

Comparison of experimental lateral separations with MC results for different zenith angle 
intervals (depth = 3750 — > 4150 m.w.e.). (D) is the average separation and Rq is the 
separation at which the decoherence function reaches its maximum. Both quantities are 
in unit of meters. 
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Table 2 

Comparison of experimental lateral separation with MC results for different rock depth 
intervals (cos^ = 0.8 — > 0.9). (D) is the average separation and Rq is the separation 
at which the decoherence function reaches its maximum. Both quantities are in unit of 
meters. 
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which introduce displacements of the muons from their original direction and position. 
The aim of this analysis is to attempt to disentangle these effects. 

For a real experiment of finite size and finite live time, it is generally impossible to 
reconstruct the full distribution, dN j dxdcf). Grillo and Parlati [16] have suggested the use 
of its first moment, 



dxd(ff ' 



which they have named as the "decorrelation function." An analytic expression for {4>{x)) 
has been derived [16] using some approximations. 

In this analysis, we have used only double muon events and employed track reconstruc- 
tion only in one projective view, namely the wire view. The results are still preliminary. 
Since the position of the shower axis is not known, the measure is relative; that is, the 
quantities (j) and x are defined between pairs of muons in the same bundle. The distance 
is always taken to be positive, while the angle is positive if the tracks are diverging, and 
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negative otherwise. Since the average angle as a function of distance is normalized to 
the number of events at that distance, the influence of apparatus effects (efficiency, work- 
ing conditions, containment) should not be important. We have made no run selection 
on the data sample, but we have made an event selection rejecting too short tracks, to 
avoid contamination of the sample by locally produced pions and small showers. The 
experimentally measured decorrelation function is presented in Figure 6. 

Figure 6 also shows the decorrelation function computed from the HFMAS code. It is 
evident that there is a strong disagreement between the data and Monte Carlo at relatively 
small distances. We have investigated several possible causes for this disagreement. Dif- 
ferent composition models have essentially no effect, nor do more refined muon transport 
codes which include Moliere tails beyond the Gaussian approximation for multiple scat- 
tering. It is possible to modify the average interaction cross section (and hence primary 
interaction height) to make the Monte Carlo results agree with the data, but the required 
modification is rather extreme and is likely inconsistent with reasonable extrapolations of 
accelerator data. We are presently investigating more subtle effects, both derived from 
the finite space resolution of the apparatus and from high energy interactions of muons 
in the rock overburden. Preliminary results are encouraging, but a more refined analysis 
is not yet complete. 
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6. CONCLUSIONS 

The large amount of underground muon bundle data collected by the MACRO experi- 
ment offer significant capability in the studies of primary cosmic ray composition and of 
interactions at very high energies. Our data do not favor the hypothesis of a dramatic 
change of primary composition towards pure Fe element immediately above the "knee." 
The discrepancy between the experimental muon pair lateral separations and the Monte 
Carlo results is under investigation. A new analysis based on the muon decorrelation 
function is being pursued. 

We are very grateful to the U.S. Department of Energy, the Laboratori Nazionali del 
Gran Sasso, and the Istituto Nazionale di Fisica Nucleare for their continued support of 
the MACRO experiment. 
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